Wigner-Crystal Formulation of Strong-Coupling Theory for Counter-ions Near Planar 
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We present a new analytical approach to the strong electrostatic coupling regime (SC), that can 
be achieved equivalently at low temperatures, high charges, low dielectric permittivity etc. Two 
geometries are analyzed in detail: one charged wall first, and then, two parallel walls at small 
distances, that can be likely or oppositely charged. In all cases, one type of mobile counter-ions 
only is present, and ensures electroneutrality (salt free case). The method is based on a systematic 
expansion around the ground state formed by the two-dimensional Wigner crystal(s) of counter-ions 
at the plate(s). The leading SC order stems from a single-particle theory, and coincides with the 
virial SC approach that has been much studied in the last 10 years. The first correction has the 
functional form of the virial SC prediction, but the prefactor is different. The present theory is 
free of divergences and the obtained results, both for symmetrically and asymmetrically charged 
plates, are in excellent agreement with available data of Monte-Carlo simulations under strong 
and intermediate Coulombic couplings. All results obtained represent relevant improvements over 
the virial SC estimates. The present SC theory starting from the Wigner crystal and therefore 
coined Wigner SC, sheds light on anomalous phenomena like the counter-ion mediated like-charge 
attraction, and the opposite-charge repulsion. 

PACS numbers: 82.70.-y, 61.20.Qg, 82.45.-h 



I. INTRODUCTION 

Understanding effective equilibrium interactions be- 
tween two charged mesoscopic bodies immersed in a so- 
lution, is essential in various fields of colloid science, 
from physics [l| to biochemistry References 
offer a general overview. A breakthrough in the field 
was achieved when it was realized in the 1980s, from 
numerical evidences, that equivalently charged surfaces 
may effectively attract each other, under strong enough 
Coulombic couplings. Such couplings can be realized in 
practic e by increasing the valency of the counter-ions in- 
volved This "anomalous" like-charge attraction ex- 
plains the formation of DNA condensates [ll| or aggre- 
gates of colloidal particles p^]. A complementary in- 
teresting although simpler to rationalize problem is the 
possibility of an effective repulsion between two plates 
with opposite uniform surface charges. 

The weak-coupling limit is described by the Poisson- 
Boltzmann (PB) mean-field approach. Formulating the 
Coulomb problem as a field theory, the PB equation can 
be viewed as the first-order term of a systematic expan- 
sion in loops U3\ . While the like-charge attraction is not 
predicted by the PB theory [l3 - [T7| . the opposite-charge 
repulsion can occur already in the mean-field treatment 
[isl [l9| , since it is merely an entropic effect with a large 
cost for confining particles in a small volume. 

A remarkable theoretical progress has been made dur- 
ing the past decade in the opposite strong-coupling (SC) 
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limit, formulated initially for a single wall or two parallel 
walls at small separation. The topic was pioneered by 
Rouzina and Bloomfield ^oj and developed further by 
Shklovskii, Levin with collaborators @, [2l| . An essential 
aspect is that counter-ions form two-dimensional (2D) 
highly correlated layers at charged walls at temperature 
T ~ 0. For small but non vanishing temperatures, the 
structure of interfacial counter-ions remains close to its 
ground-state counterpart. 

Within the field-theoretical formulation, which has 
been put forward by Netz and collaborators in [2^ , 
the leading SC behavior is a single-particle theory in the 
potential of the charged wall(s). Next correction orders 
as obtained as a virial or fugacity expansion in inverse 
powers of the coupling constant 5, defined below; we re- 
fer to this approach as the virial strong-coupling (VSC) 
theory. The method requires a renormalization of in- 
frared divergences via the electroneutrality condition. A 
comparison with Monte Carlo (MC) simulations f25] indi- 
cated the adequacy of the VSC approach to capture the 
leading large-S behaviour of the density profile, which 
was an important achievement in the field. The first cor- 
rection has the right functional form in space but an in- 
correct prefactor, whose values even depart further from 
the MC ones as the coupling constant S grows. This de- 
ficiency was attributed by the authors to the existence 
of an infinite sequence of higher-order logarithmic terms 
in the fugacity which have to be resummed to recover 
the correct value of the prefactor. The leading order of 
the VSC theory was generalized to non-symmetrically 
charges plates [l9l . f24| , image charge effects 25j , presence 
of salt [26l | and to various curved (spherical and cylindri- 
cal) geometries, for a review, see [23. Beyond Refs. [2^ . 
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several investigations assessed numerically the adequacy 
of the leading order VSC approach jH, [Ij-Iil, H, [2| . 

Since the coupling constant S cx l/T^, the zero tem- 
perature is contained in the VSC approach as the limit 
S — > oo. This question requires some care though, since 
a natural rescaled distance z = in the direction per- 
pendicular to the plate(s) is set by the Gouy-Chapman 
length /i (X T, which tends to zero as T — >■ 0. From this 
point of view, the VSC method can be seen as a low- 
temperature theory approaching T — under a special 
spatial scaling of particle coordinates. One of the restric- 
tions is the applicability of the theory to small (rescaled) 
distances between the charged plates. There exist other 
possibilities to approach the zero temperature limit. One 
of them is to construct an expansion in S around the 
limit 5 — >■ (X) , under the fixed ratio of the distance d (in 
the two plate problem) and the lattice spacing a of the 
Wigner crystal formed at T = 0. The low-temperature 
theory proposed by Lau et al. [s^l, can be considered 
in some respect as being of this kind. The considered 
model consists of two staggered hexagonal Wigner crys- 
tals of counter-ions condensed on the plates; the particles 
are not allowed to move in the slab between the plates. 
The attraction between the plates at zero and non-zero 
temperatures, which results from the interaction of the 
staggered Wigner crystals and from the particle fluctua- 
tions, can be computed. Since the particles are not al- 
lowed to leave their Wigner plane, the counter-ion profile 
between the two plates is trivial and there is no need for 
a spatial scaling. Such a model is interesting on its own, 
but has a restricted applicability to realistic systems of 
counterions because the particles are assumed to stick 
to the plates. This assumption may be perhaps accept- 
able at large distances between plates, but discards from 
the outset the excitations that are relevant at small dis- 
tances, where the counter-ions unbind from the interfaces 
(see e.g. [l^,!!^! and the analysis below). 

An interpolation between the Poisson-Boltzmann (low 
S) and SC regimes (high S), based on the idea of a "corre- 
lation hole", was the subject of a series of works (Sll - lssj . 
The correlation hole was specified empirically in Refs. 
[3^ and self-consistently, as an optimization condition 
for the grand partition function, in (ssj . An interesting 
observation in [33|], corroborated by a comparison with 
the MC simulations, was that the first correction in the 
SC expansion is proportional to 1/ a/S, and not to 1 /S 
as suggested by the VSC theory. Our exact expansion 
below shows that indeed, the first correction scales like 

Recently, for the geometries of one plate and two equiv- 
alently charged plates with counter-ions only, we pro- 
posed another type of SC theory jsj]. It is based on a 
low-temperature expansion in particle deviations around 
the ground state formed by the 2D Wigner crystal of 
counter- ions at the plate(s). The approach points to the 
primary importance of the structure of the ground state, 
a point emphasized by some authors, see e.g. l35|. Our 
starting point therefore resembles that of Ref. |30j , but in 



the subsequent analysis, the particles vibrations around 
their Wigner lattice positions are allowed along all direc- 
tions, including the direction perpendicular to the crystal 
plane along which the particle density varies in a non- 
trivial way. The theory is formulated in the set-up of 
the original VSC approach: An SC expansion around 
the same limit S — cx) is made, together with the same 
scaling of the coordinate in the direction perpendicular 
to the plate(s), z = z/fi. Since the formation of the 
Wigner crystal is the basic ingredient from which the 
method starts, we shall refer to it as the WSC theory. 
Its leading order stems from a single-particle theory, and 
is identical to the leading order obtained in the VSC ap- 
proach. In the present planar geometry, both WSC and 
VSC differ beyond the leading order, when the first cor- 
rection is considered. In this respect, in assessing the 
physical relevance of WSC and VSC, comparison to "ex- 
act" numerical data is essential. Remarkably, the first 
WSC correction has the functional form in space of the 
VSC prediction, but the prefactor is different: Its 
dependence on the coupling parameter and the value of 
the corresponding prefactor are in excellent agreement 
with available data of MC simulations, while the VSC 
prediction is off by several orders of magnitude under 
strong Coulombic couplings Unlike the VSC the- 

ory, the WSC expansion is free of divergences, without 
any need for a renormalization of parameters. The WSC 
expansion turns out to be in inverse powers of -y/S, and 
not of S like in the case of the VSC expansion. Due to 
its relatively simple derivation and algebraic structure, 
the WSC method has a potential applicability to a large 
variety of SC phenomena. In particular, the WSC can 
be worked out beyond the leading order for asymmetric 
plates, which, to our knowledge, was not done at the VSC 
level, possibly due to the technical difficulty to overcome. 
The specific 2D Coulomb systems with logarithmic pair 
interactions were treated at WSC level in Ref. fH]. 

In this paper, we aim at laying solid grounds for the 
WSC method. We develop the mathematical formal- 
ism initiated in Ref. [s^ ], which is based on a cumu- 
lant expansion, to capture systematically vibrations of 
counter-ions around their Wigner crystal positions. This 
formalism enables us to deal, in the leading order plus 
the first correction, also with asymmetric, likely or op- 
positely charged plates. The results obtained are in re- 
markable agreement with MC data, for large as well as 
intermediate values of the coupling parameter S. 

The paper is organized as follows. The one-plate geom- 
etry is studied in Sec. |lTl An analysis is made of counter- 
ions vibrations around their ground-state positions in the 
Wigner crystal, along both transversal and longitudinal 
directions with respect to the plate surface. The cumu- 
lant technique, providing us with the WSC expansions 
of the particle density profile in powers of is ex- 

plained in detail. Sec. IIIII deals with the geometry of 
two parallel plates at small separation. The cumulant 
technique is first implemented for equivalently charged 
plates and afterwards for asymmetrically charged plates. 
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In the case of the opposite-charged plates, the WSC re- 
sults for the pressure are in agreement with MC simu- 
lations for small plate separations and lead to the cor- 
rect (nonzero) large-distance asymptotics. In the case 
of the like-charged plates, the accurate WSC results for 
the pressure are limited to small plate separations. All 
obtained results represent an essential improvement over 
the VSC estimates. Concluding remarks are given in Sec. 

113 

Before we embark on our study, a semantic point is 
in order. Some authors refer to the VSC approach as 
the "SC theory" . Clearly, the VSC route is not the only 
theory that can be put forward to describe the strong 
coupling regime. In what follows, the SC limit refers to 
S — >■ cxD, and we carefully discriminate between VSC and 
WSC predictions, that will both be tested against Monte 
Carlo data. 



II. ONE-PLATE GEOMETRY 

A. Definitions and notations 

We start with the one plate problem in the 3D Eu- 
clidean space of points r = {x, y, z) pictured in Fig. [T^. 
In the half-space A' = {r, z < 0}, there is a hard wall 
of dielectric constant e which is impenetrable to parti- 
cles. A uniform surface-charge density ue, e being the 
elementary charge and cr > 0, is fixed at the wall surface 
S localized at z = 0. The g-valent counter-ions (clas- 
sical point-like particles) of charge — ge, immersed in a 
solution of dielectric constant e, are confined to the com- 
plementary half-space A = {r,2: > 0}. In this work, 
we consider the homogeneous dielectric case only, with- 
out electrostatic image forces. The system is in thermal 
equilibrium at the inverse temperature /? = \j{k^T^. 
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FIG. 1. The two geometries considered: a) one plate; b) two 
parallel plates at distance d. The neutralizing counter-ions 
have charge —qe. 



The potential energy of an isolated counter-ion at dis- 
tance z from the wall is, up to an irrelevant constant. 



E(z) 



2Trqe'^a 



(1) 



The system as a whole is electroneutral; denoting the 
(infinite) number of counter-ions by A'^ and the (infinite) 
area of the wall surface by |S|, the electroneutrality con- 
dition reads 

qN^ain. (2) 

There are two relevant length scales describing, in 
Gaussian units, the interaction of counter-ions with each 
other and with the charged surface. The Bjerrum length 

^ '-^ (3) 

e 

is the distance at which two unit charges interact with 
thermal energy /cbT. The Gouy-Chapman length 

f^=T^ (4) 

is the distance from the charged wall at which an isolated 
counter-ion has potential energy ([!]) equal to thermal en- 
ergy ksT. The z coordinate of particles will be usually 
expressed in units of fj,, 



_ z 



(5) 



The dimensionless coupling parameter S, quantifying the 
strength of electrostatic correlations, is defined as the 
ratio 



(6) 



The strong-coupling regime S ^ 1 corresponds to either 
low temperatures, or large valency q or surface charge 
ae. 

The counter-ion averaged density profile p{z) depends 
on the distance z from the wall. It will be considered in 
the rescaled form 



(7) 



27r£BCr2 ■ 

The electroneutrality condition ([2|) then takes two equiv- 
alent expressions 



q I dzp{z) 
'o 



dzp{z) 



1. 



(8) 



The contact-value theorem for planar wall surfaces [37[ 
relates the total contact density of particles to the surface 
charge density on the wall and the bulk pressure of the 
fluid P. For 3D systems of identical particles, it reads 

/3P = p{0) - 27r£Bcr^- (9) 

Since in the present case of a single isolated double layer, 
the pressure vanishes, 

2 



p{0) = 27r£BO- 



m 



1, 



(10) 



that can be viewed as a constraint that any reasonable 
theory should fulfill. 
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B. The Virial Strong Coupling approach 

With our choice of reduced units, the exact density 
profile is a function of two variables only: S). It 
is well behaved when S — >■ oo, which is nevertheless a 
limit where in unsealed variables, all counterions stick 
to the plate, forming the Wigner crystal {p{z, S) oc 5[z) 
for S — >• oo). The purpose of the present discussion is to 
resolve the structure of the double-layer at large but finite 
S. According to the VSC method 0, [13, the density 
profile of counter-ions can be formally expanded in the 
SC regime as a power series in 1/S: 



100 r 



p(z,^) = po(2) + l:Pl{z 



where 



Po{z) 



(11) 



(12) 



The leading term po(z), which comes from the single- 
particle picture of counter-ions in the linear surface- 
charge potential, is in agreement with the MC simula- 
tions [22]. Indeed, for large S, the particles' excursion 
perpendicular to the plane, which is always quantified 
by /X, is much smaller than the lateral spacing between 
ions (denoted a below) As a consequence, these ions 
experience the potential of the bare plate, while the inter- 
actions with other ions become negligible by symmetry. 
On the other hand, the MC simulations indicate that the 
sub- leading term pi (z) has the expected functional form 
(for sufficiently large coupling S > 10), but the prefactor 
1/S is incorrect. On the basis of the prediction ([TH) . the 
MC data were fitted in [22| by using the formula 



Piz,^) - Po{z) = -pi(F), 



(13) 



where p(z, S) is the density profile obtained from MC 
simulations and 9 is treated as a fitting parameter. Ac- 
cording to the VSC result (fTTj) . 9 should be given by 
= S plus next-leading corrections. As is seen in the 
log-log plot of Fig. [21 the numerically obtained values of 
9 are much smaller than S, and the difference between 9 
and S even grows with increasing the coupling constant. 



C. The Wigner Strong Coupling expansion 

Our approach is based on the fact that in the asymp- 
totic ground-state limit 5 — cxd, all counter- ions collapse 
on the charged surface z = 0, forming a 2D Wigner 
crystal ^ [2l|- It is well known that the lowest 
ground-state energy for the 2D Wigner crystal is pro- 
vided by the hexagonal (equilateral triangular) lattice. 
Each point of this lattice has 6 nearest neighbors form- 
ing a hexagon, see Fig. [3] The 2D lattice points are 
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FIG. 2. The fitting parameter 6, defined by Eq. (fT3)) . vs. the 
coupling constant H for one-plate geometry. The MC values 
reported in Ref. [S^] are shown with filled diamonds, the 
original prediction S = E of the VSC theory with the dashed 
line; the sohd curve is for our WSC prediction, given by Eq. 



indexed by {j = (ji, ^'2)}, where ji and j2 are any two 
integers (positive, negative or zero): 



where 



= iR^,Rj)=jiai+j2a2, 



/ 1 Vs' 

ai=a(l,0), "2 = a|2'~ 



(14) 



(15) 



are the primitive translation vectors of the Bravais lattice 
and a is the lattice spacing. Since at each vertex, there 
is just one particle, we can identify j with particle labels, 
j = 1,...,N {N — >■ cxi). There are two triangles per 
vertex, so the condition of global electroneutrality ^ 
requires that 



q 

a 



n/3 



(16) 



Note that in the large-S limit, the lateral distance be- 
tween the nearest-neighbor counter-ions in the Wigner 
crystal a is much larger than the characteristic length p 
in the perpendicular z-direction, aj p oc -x/S ^1. As 
invoked above, this very feature explains why a single 
particle picture provides the leading order term in a SC 
expansion, so that the two different approaches discussed 
here (VSC and WSC) coincide to leading order. The 
same remark holds for the two plates problem that will be 
addressed in section IIIII It should be emphasized though 
that this coincidence of leading orders is specific to the 
planar geometry. The z-coordinate of each particle in the 
ground state is zero, Zj = 0. 

We denote the ground-state energy of the counter- 
ions on the Wigner lattice plus the homogeneous surface- 
charge density ue by Eq. For S large but not infinite, the 
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FIG. 3. Hexagonal structure of the 2D Wigner crystal: ai 
and 02 are the primitive translation vectors. 



fluctuations of ions around their lattice positions, in all 
three spatial directions, begin to play a role. Let us first 
shift one of the particles, say j — 1, from its Wigner lat- 
tice position (Ri, Zi = 0) by a small vector Sr — {x, y, z) 
{\5r\ ^ a) and look for the corresponding change in the 
total energy SE = E — Eq > 0. The first contribution to 
SE comes from the interaction of the shifted counter-ion 
with the potential induced by the homogeneous surface 
charge density: 



SE^^\z) 



2'iTqe'^a 



(17) 



The second contribution to 5E comes from the interac- 
tion of the shifted particle 1 with all other particles j ^ 1 
on the 2D hexagonal lattice: 



5E^^\x,y,z)^ 
[qef 



E 



+ (^1, + y? + 



1 



(18) 



where = {R^^,Ry) = Ri - Kj and Rij = |R 



Rescaling the lattice positions by a and taking into ac- 
count the inequalities x/a,y/a, z/a <^ 1, this expression 
can be expanded as an infinite series in powers of x/a, 
y/a and z/a by using the formula 



1 



VTTt 



1 



1 



-t + -r 
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t < 1. 



Up to harmonic terms, the expansion reads 

'1 



5i?(^)(x,y,z) = M!c3 



-{x +y )-z 



(19) 



(20) 



Here, C3 is the special s — 3 case of dimensionless hexag- 
onal lattice sums 



a. 



1 



which can be expressed from the general theory (39| as 



E 



3,fc=-oo 
(3.fc)75(0,0) 



iP + jk + fc2)3/2 



7=3- 



c 



(22) 



with C(z, q) — J2'^=o 1/(9 + "•)^ the generalized Riemann 
zeta function and (^(z) = (^{z, 1) (this function should 
not be confused with the parameter C, appearing with- 
out arguments below after Eq. ((63)) . that will measure 
the asymmetry between two charged plates). Explicitly, 
C3 = 11.034.... The absence of the linear x, y terms 
and of the mixed xy term in (|20p is caused by the fact 
that every lattice point is at a center of inversion. The 
invariance of the hexagonal lattice with respect to the 
rotation around any point by the angle 7r/3 implies the 
lattice sum equalities 



1 

2 



9 E /(^ii)-^ij ' 



(23) 



which were also used in the derivation of pO|. Note 
that the x'^ and j/^ harmonic terms in Eq. ([20]) have 
positive signs which is consistent with the stability of 
the Wigner crystal in the {x, y) plane. On the other 
hand, the minus sign of the 2^ term does not represent 
any stability problem due to the presence of the posi- 
tive linear contribution in (|17|) . which is dominant for 
small z-distances. The total energy change is given by 
SE{x,y,z) = SE^^^z) + SE^^^x,y, z). Finally, let us 
write down the z-dependent part of the dimensionless 
energy shift —PSE, with z expressed in units of /i: 



|35E{Q,Q,^xI) - -z + 



3V4 
2^' 



(24) 



We see that in the limit S — > 00, as advocated above, 
the two-body interaction term of the shifted ion with all 
other ions on the Wigner crystal is of order l/-\/S and 
therefore negligible in comparison with the one-body po- 
tential term —z due to the surface charge density. This 
leading single-particle picture is common to both VSC 
and WSC approaches. As concerns the two-body inter- 
action terms of higher orders (p = 3, 4, . . .), their co- 
efficients are proportional to (f' t^\x^ j dP^^ cx I/S'p^^^/^. 
The present scheme thus represents a systematic basis 
for an expansion in powers of 1/a/S. 

The generalization of the above formalism to inde- 
pendent shifts of all particles from their lattice posi- 
tions is straightforward. Let us shift every particle 
j ~ 1,2,...,N from its lattice position (ILj,Zj — 0) 
(21) by a small vector Svj = {xj,yj,Zj) {\Srj \ <C a) and study 
the corresponding energy change SE. As before, the first 
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(one-body) contribution to SE is given by 



N 



i=i 



(25) 



The second (two-body) contribution to SE is expressible 



as 



6E^^H{x,}Ayj}A^j}) 



N 



3,fc = l 



^1 + fijk + Vjk 



(26) 



where the dimensionless jij^ and Ujk involve the parti- 
cle coordinates along and perpendicular to the Wigner 
crystal, respectively: 



TDX 

Hjk = 2{xj -~Xk)-^ + 2{yj - Vk) 2 



[i^i ~ ^kf + {yj - Vkf] , (27) 



■jk 



1 



= -^(Zj - Zk) 



(28) 



Performing the expansion of type (|19p in small /ijfc and 
Vjk^ we end up with 

^ pSE^^^ixAAvAAzj}) ^ + Sw + S,,w, (29) 



where 



(30) 



contains particle shifts exclusively in the z direction, 

N 



Sw — 



E 



1 /I 



2 iJjfe V2^^^ 8^-'* 



j.k = l 



(31) 



contains particle shifts exclusively in the (a:, y) Wigner 
plane and 



5, 



z,W 



q^h 
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^ 1 

F — 



+Y6 il^h"^'' + t^^'^'^Jk) + 



(32) 



mixes particle shifts along the z direction with those 
along the (x, y) plane. 

We are interested in the particle density profile de- 
fined by p(r) = (X^jLi '^(j' ~ where (• • •) means 
thermal equilibrium average over the Boltzmann weight 
exp{—(3SE) with 



/3SE = ~I3SE^^^ ~ PSE^^^ 

N 

= — ^ Zj + Sz + Sw + Sz^w- 

3 = 1 



(33) 



The ground-state energy Eq is a quantity which is in- 
dependent of the particle coordinate shifts and as such 
disappears for the statistical averages. The system is 
translationally invariant in the (a;, y) plane, so that the 
particle density is only z-dependent, p{r) — p{z). We 
shall consider separately in ([33]) the terms containing ex- 
clusively particle shifts in z direction, transversal to the 
wall, and those which involve longitudinal particle shifts 
along the Wigner [x, y) plane. 



D. Contribution of transversal particle shifts 

Let us forget for a while the terms Sw and Sz,w in 
([33)1 and consider only the particle z-shifts in the "most 
relevant" Sz, 



(36E 



N 

-Y.z,+Sz 



(34) 



Expressing z in units of /i, 5*2 in Eq. (|30p can be written 
as an infinite series in powers of l/\/S, the first terms of 
which read 



Sz = 



a 



N 

E 



1 



4VS {Rjk/a) 

(j^k) 

N 



3a^ 



E 



16S3/2 {R,k/a) 



3 (Zj - Zk) 

(zj - Zkf + 



(35) 



In the limit S ^ oo, 5*2 is a perturbation with respect to 
the one-body part in (p4| . 

To obtain the particle density, we add to the one-body 
potential z an auxiliary (generating or source) potential 
/3u(r), which will be set to at the end of calculations. 
The partition function of our A'^-particle system 



1 I- ^ 



driw(ri)e exp{Sz) (36) 



thereby becomes a functional of the generating Boltz- 
mann weight w{r) = exp[— /3u(r)]. The particle density 
at point r is obtained as the functional derivative 



P(r) 



6w{r) 



In Zpf [w] 



w{r) — l 



(37) 



which is of course a function of S, in addition to r. To 
treat Sz as the perturbation, we define the Sz — coun- 
terpart of the partition function p6p 



drw{r)e~ 



1 

m 



(38) 
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which corresponds to non-interacting particles in an ex- 
ternal potential. It is clear that 




= ln(exp(5^))o, 



(39) 



where (• • •)o denotes the averaging over the system of 
non-interacting particles defined by We are left 

with the cumulant expansion of ln(exp(>5'3))o: 

oo 

1 



n=l 



.(40) 



An important property of the cumulant expansion is that 
if {Sz)o is an extensive (proportional to N) quantity, the 
higher-order terms will also be. In other words, the con- 
tributions of N'^, N^, etc. orders will cancel with each 
other. We conclude that 

\nZN[w] = InZ^^^^H + (5,)o -t- ^ {{Sl)o - {S,)l) +■■■. 

(41) 

The particle density results from the substitution of this 
expansion into p7p. and the subsequent application of 
the functional derivative with respect to w(r), taken at 
w{r) = 1. 

The leading SC behavior of the particle density stems 
from In [w] . Since 



Sw{r) 



w{r)=i Jj^dre-^ " \j:\fi 
= (27r^BO-2)e-" 



(42) 

we have po{z) ~ e~^, which coincides with the leading 
VSC term presented in (|12|) . 

The first correction to the density profile stems from 
{Sz)o, namely from the first term in the series represen- 
tation of Sz 



{Sz)o 



a' 



N 



1 



4VS f^, (Rjk/a) 



3((2|+2^-2J,i-,))o. (43) 



It is easy to show that 
6 



6w{r) 



w{r) — l 



-z (jp 



(46) 



where we used the equality [2'']oU(r)=i — P- The for- 
mula for the density profile, in the leading order plus the 
first correction, then reads 



p(z,s) 



3^/4 Cs 



z]+0 



(47) 



Note that the electroneutrality ([SJ and the contact the- 
orem pO| are satisfied by this density profile. In Fig. |4l 
we compare the appropriately rescaled first correction to 
the leading SC profile obtained in (|T7)) (solid curve) with 
MC data 22[ at S = 10"^ (filled squares). The agreement 
is excellent. On the other hand, the VSC prediction is 
off by a factor lOOO^/^. 



-0.1 




FIG. 4. Single charged wall: Comparison between the rescaled 
analytical first correction to the strong coupling profile from 
Eq. gZl) (solid curve) and the MC results of Ref. ^ (filled 
squares). Here, H = 10^ and po{z) denotes the leading order 
term exp(— i), that is subtracted from the numerical data to 
probe the correction. 



A useful property of the averaging (• • ■)o is its inde- 
pendence on the particle (lattice site) index, e.g. for 
p = 1, 2, ... we have 

IaUZi [dr^w{r,)e 



IAIl^=l [dr^w{r,)e 
J^dTw{r)e- 



■z~p 



drw{r)e^^ 



(44) 



Simultaneously, due to the absence of interactions in 
(• • •)o, correlation functions of particles decouple them- 
selves, e.g. {zjZk)o — [2]o for j ^ k. Thus, the relation 
becomes 



2VS 



=7VC3 ([i^]o - [i] 



(45) 



Comparing our WSC resuft (gT]) with the VSC Eqs. 
(fTTjl and (fT2|) we see that the first corrections have the 
same functional dependence in 2, but different pref actors. 
In terms of the fitting parameter 9 introduced in (|13p . the 
VSC estimate = S is compared with the present value 



8^3/2 1 
33/4 C3 



S = 1.771, 



(48) 



As is seen from Fig. [21 this formula (solid curve) is in 
full agreement with the data of MC simulations (filled 
diamonds) . 

In the series representation of Sz (|35p . the first term is 
of order S"^/^ and the second one is of order S~3/2 
view of (|¥T|) . the second correction to the density profile 
stems from ((S'f)o — {Sz)'q)I'2- with Sz represented by its 
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first term, and not from {Sz)o with Sz represented by its 
second term. Let us analyze in detail the average 



\ y-^ 1 X - 1 



(49) 



For a fixed pair of site indices (j ^ A;), there exist seven 
topologically different possibilities for the pair (to ^ 7i): 



TO = J, 71 = fc; 
n = j, TO = k\ 



n = j, TO 7^ j, /c; 

m = k, n ^ j, k; 

n = k, m ^ j, k; 

mj^j,k, j, k, to.} 



factor 2 



factor 4 



factor 1 



Here, respecting the properties of the averaging (• • •)o, 
those possibilities which lead to the same result are 
grouped together. After simple algebra, we find that 



{S'z)o - ^{nc! - 4[z%[2]o + 3[zro) 



ANCi + 2NCq 



(50) 



The "undesirable" disconnected term of order N'^ is can- 
celled by the subtraction of {Sz)o- After performing 
the functional derivatives with respect to u'(r), taken at 
^(r) — 1, we end up with the next correction to the 
profile of the form 



33/2 1 

647r3 

(51) 

Note that this correction does not break the electroneu- 
trality condition (jS)) nor the contact theorem pO)) . 



E. Contribution of longitudinal and mixed particle 

shifts 



Now we consider in ([33l) also the term Sw with purely 
longitudinal particle shifts in the Wigner plane and the 
term Sz,w with mixed transversal and longitudinal shifts. 
Denoting particle shifts in the infinite Wigner plane as 
Uj = {xj,yj), these terms possess the important transla- 
tional symmetry: 



Sw{{n,}) = Sw{{uj+u}), 

Sz,w{{Uj,Zj}) = Sz.wii^] + U> Zj}), 



(52) 



where u is any 2D vector. We first investigate the scaling 
properties of Sw and Sz.w- 

Let us expand Sw up to quadratic x, y-deviations: 

^ _ q'iB ^ {R%/a)'-2{R^Jar f 



N 

E 



{R^^Jaf-2{R%laY ( - y. 



4a {Rjk/ar 



3q2 



N 

E 



{R%R%)l0''^ {xj - xk){yj ~ Vk) 



2a {Rjk/af 



+ ■ 



(53) 



Terms linear in {xj — Xk)/a and {yj—yk)/^ vanish because 
every point of the hexagonal Wigner crystal is a center 
of inversion. We saw that in the z direction the rele- 
vant length scale is determined by the Gouy-Chapman 
length /i: Rescaling the z coordinates by /i, the (leading) 
linear potential term z is independent of the coupling 
constant S while the next terms are proportional to in- 
verse powers of \/S and therefore vanish in the SC limit. 
The natural length scale in the Wigner {x,y) plane is 
the lattice spacing a, but this is not the relevant scale 
in statistical averages. The relevant length A is deter- 
mined by the requirement that the rescaling of coordi- 

XYj in ([55)) leads to a di- 



nates Xj 



XXj and yj 



mensionless and independent (leading) quadratic term. 
Since g^fe/a « we have 



A 1 

— C>C ^777. 



A 



(54) 



(the numerical prefactors are unimportant), i.e. the rel- 
evant scale is "in between" fi and a. The higher-order 
terms in Sw, which contain the deviations {xj — x^) and 
iVj ~ Vk) in powers p = 3, 4, . . ., scale like l/S^^*"^'/* and 
therefore vanish in the limit S — cx). 

Let us now consider the leading expansion terms of the 
mixed quantity Sz,w- 



S: 



+ 



Sq^iB 


N 

E 


[{zj - Zk)/af 


4a 


{Rjk/af 






-^k\^ R% 


a 




a J a 


Sq^iB 


N 

E 

j,k = l 


[{zj - Zk)laf 


8a 


{Rjk/ay 



V] - Vk 



R 




Xk 




R 



Vj - Vk 
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+10- 



a 



a 



+ •••. (55) 

Rescaling the particle coordinates as follows Zj = fiZj, 
Xj = XX j , yj = XYj , the first term is of order 1 ja?!'^ and 
the second one is of order 1/S. 

To obtain the density profile, one proceeds in analogy 
with the previous case of transversal vibrations. We in- 
troduce the partition function of our A^-particle system 



1 

iV! 



N 



(56) 



with the generating Boltzmann weight w(r). We take as 
the unperturbed system the one with one-body potentials 
— Zi in z direction and Sw in (x, y) plane, and treat 5*^ + 
Sz,w as the perturbation. Using the cumulant method, 
we obtain 



lnZAr[w] = \wZ^^\w\ + (S'^)o + (S'z,vy)o + 



(57) 



where (• • •)o denotes the averaging over the unperturbed 
system with the partition function 



1 r ^ 



dviw(j:i)e ^' exp(S'vK)- (58) 



The particle density is given by Eq. ([57]) . 

The additional appearance of exp(S'vi/) in the averaging 
over the unperturbed system is a complication which can 
be sometimes removed trivially by using the translational 
invariance of Sw (|52p . We shall document this fact on 
the leading SC behavior of the particle density at point 
r — (u, z) which stems from lnZ^''[w]: 



'w{y') — \ 



S-w 



(59) 



Since the surface of the plate S is infinite, we shift in 
the denominator the integral variables z 7^ 1 as follows 
Ui — U; -t- Ui — u which transforms Sw — Sw(\\\ = u). 
Integrating over Ui, the ratio of integrals in ([5^ is u- 
independent, and reads By this simple technique, 

it can be shown that the contribution to the density pro- 
file coming from the functional derivative of (52)0 is not 
affected by Sy/-, which decouples from the z- variables. 
We remember from the previous part about transversal 
deviations that (S'z)o is of order 

The description is a bit more complicated in the case 

of 



/Alli^i [^i"i'^(i"i)e exp( 5*1^)5; 



/a Dili Mr,w(r,)e-^.]exp(S', 



(60) 



In the corresponding contribution to the density profile, 
obtained as the functional derivative with respect to w(r) 
at w{t) — 1, the z and (x, y) subspaces decouple from one 
another. The z variables are considered in the rescaled 
form z — z/pL. To perform the integration over the 
Wigner plane, we rescale the (a;,y) variables to the ones 
\{X, y); this ensures that the quadratic part of Sw is 5- 
independent and all higher-order terms p = 3, 4, . . ., pro- 
portional to l/S^^*"^)/**, vanish in the SC limit S 00. 
Thus the leading dependence on 2 is given by the scal- 
ing factor of Sz^w under the coordinate transformations 
z — fiz and {x,y) — X{X,Y), which was found to be of 
order l/S'^/^. This contribution does not alter the first 
correction cx l/\/S. To calculate explicitly the second 
correction is a complicated task, because the quadratic 
part of in the exponential eiq){Sw) involves all inter- 
actions of particles on the Wigner crystal. The explicit 
diagonalization of Sw can be done e.g. in the small wave 
vector limit [38j. 

The fact that the longitudinal vibrations in the plane 
of the Wigner crystal have no effect on the leading term 
and the first correction of the particle density profile is a 
general feature of the WSC theory. In what follows, we 
shall ignore these degrees of freedom, restricting ourselves 
to the leading term and the first correction, proportional 
to 



PARALLEL PLATES AT SMALL 
SEPARATION 



III. 



Next we study the geometry of two parallel plates 
El = 1 and E2 = 2 of the same (infinite) surface 
|Ei| — IE2I = |E|, separated by a distance d, see Fig. 
[TId. The 2 = plate 1 carries the constant surface charge 
density trie, while the other plate 2 at z = d is charged 
by (726. The electric potential between the plates is, up 
to an irrelevant constant, given by 



(^) 



27r(cri - (T2)e 



(61) 



N mobile counter-ions of charge —qe (the valency q > 
0), which are in the region between the walls A = 
{r,0 < z < d}, compensate exactly the fixed charge on 
the plates: 



giV=(ai+(72)|E|. 



(62) 



Without any loss of generality we can assume cri > 0, so 
that the asymmetry parameter 



c 



CT2 
CTl 



> 



(63) 



w) 



This parameter should not be confused with the Riemann 
function introduced in Eq. (j22[) . By rescaling appropri- 
ately model's parameters, it is sufficient to consider the 
interval — 1 < C ^ 1- The limiting value C — corre- 
sponds to the trivial case a2 = —ci with no counter- ions 
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between the plates. The symmetric case C = 1 corre- 
sponds to equivalently charged plates 0-2 = fi . Note that 
in all cases considered, there is only one type of mobile 
ion in the interstitial space < z < c?. 

Because of the asymmetry between the surface charges, 
there exist two Gouy-Chapman lengths 



1 



1 



Ail = 



(64) 



2neBqcT,-^' 27r^Bg|a2| |CI ' 

Similarly, we can define two different coupling parameters 



Ml 



M2 



IC|S. 



(65) 



Here, for the ease of comparison, we follow the convention 
of Ref. 01 : all quantities will be rescaled by their plate 
1 counterparts, i.e. z = z/jii, and 



P 



(66) 



The reduced density is a function of three arguments: z, d 
and S while the reduced pressure depends on two: d and 
S. For notational simplicity, the dependence on d and 
S will often be implicit in what follows. Note also that 
P = eP/(27re^(T^), so that the rescaling factor required 
to defined the dimensionless pressure is temperature in- 
dependent. This is not the case of the rescaling factor 
applied to distances, since the Gouy-Chapman lengths 
scale as T. The electroneutrality condition (|62l) can be 
written in two equivalent ways 



dzp{z) 



CTl + 0-2 



dzpC^ = 1 + C- (67) 



The contact-value theorem ([S]), considered at z = and 
z = d boundaries, takes two equivalent forms 



p(o) - 1 = p(d) - c^ 



(68) 



which provides a strong d and independent constraint 
for p(0) - p{d). 

In the case of oppositely charged surfaces — 1 < C ^ Oi 
the ground state of the counter-ion system is the same as 
for the isolated plate 1, i.e. all N counter- ions collapse 
on the surface, and create the hexagonal Wigner crystal. 
For this region of C values, one can easily adapt the WSC 
technique from the one-plate geometry for a priori any 
distance d between the plates. 

The case of like-charged plates < C ^ 1 is more 
subtle. The ground state of the counter-ion system cor- 
responds to a bilayer Wigner crystal, as a consequence of 
Earnshaw theorem [40] . The lattice spacings of each layer 
are denoted hi and 62; they are the direct counterpart of 
the length scale a introduced in section |lll The bilayer 
structure i s, in ge neral, complicated and depends on the 
distance d |4l| - t43| . For this region of C, values, the WSC 
technique cannot be adapted directly from the one-plate 



geometry, except for small distances between the plates 
such that d <^ b, where b — min{foi, 62}- The point is that 
each particle experiences, besides the direct linear one- 
body potential (ISTI) induced by homogeneously charged 
plates, an additional perturbation due to the repulsive 
interactions with other g-valent ions. This additional po- 
tential is, for d <^b, small compared to (pT|) . This opens 
the way to a perturbative treatment along similar lines 
as in section ini in which the leading one-body descrip- 
tion is then fully equivalent to the one derived within the 
VSC method. 

First we shall address the symmetric C = 1 case which 
ground state was studied extensively in the past. The 
symmetric configuration is of special importance in the 
VSC method: Although the leading SC result for the 
density profile and the pressure was derived for all val- 
ues of the asymmetry parameter — 1 < C ^ 1 Hi; the 
first SC correction (inconsistent with MC simulations) is 
available up to now only for C = 1 [IS 113 ■ After solv- 
ing the SC limit for the symmetric case, we shall pass to 
asymmetric, oppositely and likely charged, surfaces and 
solve the problem in the leading SC order plus the first 
correction. 



A. Equivalently charged plates 

For (71 = CT2 = (7, the electric field between the walls 
vanishes. At T = 0, the classical system is defined fur- 
thermore by the dimensionless separation 
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1 



/2^^/s' 



(69) 



A complication comes from the fact that counter-ions 
form, on the opposite surfaces, a bilayer Wigner crystal, 
the structure of which depends on rj [4l|-|43|. Two limit- 
ing cases are clear. At the smallest separation r/ = 0, a 
single hexagonal Wigner crystal is formed. Due to global 
neutrality, its lattice spacing h is given by 



-1 = 

2a 2 ' 



(70) 



The lattice spacing is simply related to that of the one 
plate problem by 6 = a/\/2. At large separations 77 — >■ cx), 
each of the plates has its own Wigner hexagonal struc- 
ture and these structures are shifted with respect to one 
another. The transition between these limiting phases 
corresponds to the following sequence of structures (in 
the order of increasing 77 |41|): a mono-layer hexago- 
nal lattice (I, < 77 < 770), a staggered rectangular 
lattice (II, 770 < 77 < 0.26), a staggered square lattice 
(III, 0.26 < 77 < 0.62), a staggered rhombic lattice (IV, 
0.62 < 77 < 0.73) and a staggered hexagonal lattice (V, 
0.73 < 77). The three "rigid" structures I, III and V, 
which do not change within their stability regions, are 
shown in Fig. [S] The primary cells of intermediate "soft" 
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II and IV lattices are changing with rj within their stabil- 
ity regions. The existence of phase I in a small, but finite 
interval of ?7, is a controversial issue [4TI - l43l [. and there- 
fore, so is the case of the precise value of the threshold 
r]Q. Whether rjQ is vanishing or is a very small number, 
remains an open problem. Here, we perform expansions 
of thermodynamic quantities in powers of d/b <^ 1 (or, 
equivalently, r/ cx d/VE < 1 since the scale d is fixed 
while S becomes large). We therefore need to know the 
ground state structure for d/b oc 77 = 0, which is clearly 
structure I, irrespective of the "770 controversy" , with a 
lattice spacing given by ((70|) . We shall thus document 
our WSC expansion on structure I. 



Structure I 



thus reads 



rp(2) i<ie? (Zj -Zkf 



+ 



d^ 



Ae 

2e IR R |3- 



(71) 



Note that the first (quadratic in z) term carries only the 
information about the single Wigner crystal of lattice 
spacing b. The information on how the lattice sites are 
distributed between the two plates within structure I is 
contained in the second constant (from the point of view 
of thermal averages irrelevant) term which compensates 
the first one if the counter-ions are in their ground-state 
configuration. The harmonic terms in the (x, y) plane 
prove immaterial for the sake of our purposes. The total 
energy change is given, as far as the z-dependent contri- 
bution is concerned, by —f36E = Sz with 



Sz 



4VS 



(72) 
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FIG. 5. Rigid ground-state structures I, III and V of counter- 
ions on two parallel charged plates; open and filled symbols 
correspond to particle positions on the opposite surfaces. 



Let Rj = (i?J,i?p be the position vector of the par- 
ticle localized on the shared hexagonal Wigner lattice of 
type I; Zj = if the particle j — 1, . . . , N/2 belongs to 
the plate Ei (say filled symbols of Structure I in Fig. 
O and Zj = if the particle j = N/2 -t- 1, be- 
longs to the plate S2 (open symbols of Structure I in 
Fig. [S|). Let us shift all particles from their lattice posi- 
tions {TLjjZj = V d} to {{xj,yj, Zj)} and look for the 
corresponding energy change SE from the ground state. 
Since the potential induced by the surface charge on the 
walls is constant between the walls and the linear in z 
contribution of Wigner crystals is negligible if d/b ^ 1, 
the corresponding SE^-^^^ — 0. The z-coordinates of par- 
ticles, constrained by the distance d between the plates, 
are much smaller than the Wigner lattice spacing b, i.e. 
both and {zj — ZkY are much smaller than |Rj — R^P 
for j ^ k. The harmonic in z part of the energy change 



The only difference between this two-plate Sz and the 
one-plate Sz ((55)) consists in the factor 2'^/^ due to the dif- 
ferent lattice spacing of the corresponding Wigner crys- 
tals, b = a/\f2. 

To derive the density profile, we use the cumulant tech- 
nique with the one-body Boltzmann factor equal to 1 (no 
external potential). The leading SC behavior stems from 



{^^dvw(r:)\^ /N\. Since 



(5w(r) 



N 

w{v)=\ \Y\d 



{2^iwy- 



(73) 



we have in the leading SC order the constant density 
po (2) ~ 2/ d. This is the one-particle result in zero poten- 
tial, respecting the electroneutrality condition (I57|) with 
C = 1. The same leading form was obtained by the VSC 
method [2^ [2^. The physical meaning is simple: due 
to their strong mutual repulsion, the counter-ions form a 
strongly modulated structure along the plate and conse- 
quently decouple in the transverse direction, where they 
only experience the electric field due to the two plates. 
In the symmetric case C = 1, this field vanishes and 
the resulting ionic density is uniform along z: from elec- 
troneutrality, it reads po — 2/d. The situation changes in 
the asymmetric case, where one can anticipate poi again 
driven by the non vanishing but uniform bare plates field, 
to be exponential in z. 

The first correction to the density profile stems from 



(5. 



-j^NCi {[z% 



where 



K 



J^drw{r) 



p= 1,2,, 



(74) 



(75) 
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Simple algebra yields 
S 



6w{r] 



1 



Ji>(r) = l 



P+l 



(76) 



where we used that [2:^]oU(r)=i = d.^/{p+^)- The density 
profile is thus obtained in the form 



2 12 




12 



, (77) 



where 



9{C = 1) 



(477)^/^ 1 1 

33/4 C3 V2' " 



S = 1.252. 



(78) 



This density profile respects the electroneutrality condi- 
tion (157)) with C = 1- The functional form of ([77)) coin- 
cides with that of Moreira and Netz 22, 23]. For (not yet 
asymptotic) S — 100, the previous VSC result 9 = E is 
far away from the MC estimate 9 ~ 11.2 [22|, while our 
formula (|78)) gives a reasonable value ~ 12.5. 

In the evaluation of the 9 factor in Eq. (|75)) . we 
use the exact result (E^ for the lattice sum C3 of the 
mono-layer hexagonal structure I, which was the start- 
ing point of our expansion. It is instructive to compare 
(|75)) with the corresponding 9 factors calculated for the 
structures III and V presented in Fig. [5l Using a rep- 
resentation of the lattice sums in terms of quickly con- 
vergent integrals over products of Jacobi theta functions, 
we find that 9 = 1.232 . . . \/S for the structure III and 
9 = 1.143 . . . for the structure V. These values show 
only a slight dependence of 9 on the structure of the 
ground state. 

Applying the contact-value theorem (|68)) to the density 
profile ([77)) . the pressure P between the plates is given 
by 



P = -1 



39 



O 



(79) 



A similar result was obtained within the approximate 
approach of Ref. 33], with the underestimated ratio 
6i/%/S = 3V3/2== 0.866.... 

Equation (|7^ provides insight into the like charge at- 
traction phenomenon. The attractive {P < 0) and repul- 
sive {P > 0) regimes are shown in Fig. [51 Although our 
results hold for d <^ and for large S, the shape of the 
phase boundary where P = (solid curve) shows strik- 
ing similarity with its counterpart obtained numerically 
22,[32l- -l^or instance, the terminal point of the attraction 
region, shown by the filled circle in Fig. |6l is located at 
d = 4, a value close to that which can be extracted from 
[2^ [33 ]. However, for S < 20, our results depart from 
the MC data, and in particular, WSC underestimates the 
value of S at the terminal point: we find Sterm — 4.53 
(corresponding to a critical value Sjerm =8/3), whereas 
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FIG. 6. Phase diagram following from the WSC equation of 
state ((79}, for symmetric like-charged plates {C ~ 1). The 
solid curve, which shows the points where P = 0, divides 
the (H, d) plane onto its attractive {P < 0) and repulsive 
(P > 0) parts. The dashed line is the original VSC prediction 
[23 |. The filled squares are the MC data from Ref. [2^ with 
E > 20. The filled circle indicates the terminal point of the 
attraction/repulsion separatrix, obtained within WSC. The 
question mark is a reminder that the upper branch of the 
isobaric curve P = is such that d oc \/H, whereas our results 
are meaningful under the proviso that d ^ ^/E,. 



the numerical data reported in [22| yields Eterm — 12. 
The previous results apply to the VSC approach as well, 
where the functional form of the equation of state is the 
same as in WSC. Since we have 6* = S in VSC, we con- 
clude that Eterm = 8/3 ~ 2.66 within VSC, which is 
indeed the value that can be seen in Fig. [S] Clearly, 
accounting correctly for the behaviour of the counter-ion 
mediated pressure for S < 20 requires to go beyond the 
strong-coupling analysis. In addition, one has to be cau- 
tious as far as the location of the upper branch of the 
attraction/repulsion boundary is concerned: It is such 

that d/^/E is of order unity and hence lies at the border 
of validity of our expansion. 

There is another feature of the equation of state un- 
der strong coupling that can be captured by our analy- 
sis: The distance of maximal attraction, where the pres- 
sure is most negative. We predict the maximum at- 
traction, following from dP/dd — 0, to be reached at 
dmax = VQ9 (X E^^'^. Since dma^/VE (X 5~^/4 _^ g in 
the asymptotic limit S — >■ cx), we can consider the latter 
prediction, shown by the dashed line in Fig. [TJ as asymp- 
totically exact. We note that it is fully corroborated by 
the scaling laws reported in [s^, while VSC yields the 
scaling behaviour dmax S^/^. 

We now analyze in more details the short distance be- 
haviour of the pressure. The difference P — 2/d, which 
is equal to —1 in the leading SC order and is linear in 
d as concerns the first correction, is plotted in Fig. [8| 
as a function of the (dimensionless) plate separation d. 
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FIG. 7. The symmetric case = 1: The maximum attrac- 
tion distance dmax (dashed Une) is defined by dP/dd — 0. 
The sohd curve d* is the boundary between attractive and 
repulsive regimes. 



ment with the MC simulations [l^. It is interesting to 
note that in the distance range d < 2, the 5 = 0.5 data 
depart from the mean-field (Poisson-Boltzmann) results 
[23 ]. see Fig. [8j there, the inter-plate distance becomes 
comparable or smaller to b, which means that the dis- 
crete nature of the particles can no longer be ignored; 
At larger distances only does the continuum mean-field 
description hold. For small inter-plate distances, we ex- 
pect the single particle picture to take over, no matter 
how small S is. This explains that P — 2/d ^ —1, but 
there is then no reason that WSC or VSC would provide 
the relevant d correction at small S. The fact that WSC 
and VSC agree with each other here at S = 0.5 is a hint 
that such a correspondence with MC is incidental (and 
indeed, in this range of couplings, S and S^/^ are of the 
same order). It would be interesting to have MC results 
at very small S values, and to concomitantly develop a 
theory for the first pressure correction to the leading term 
2/d- 1. 
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FIG. 8. The dependence of P — 2/d on the plate separation 
d for three values of the coupling constant H = 100, 10 and 
0.5. Here C = 1 (symmetric case). The plots yielded by the 
WSC equation of state (|79p are represented by dashed lines. 
Monte Carlo data are shown with symbols: open circles 
for H = 100, filled diamonds for H = 10 and open diamonds for 
H = 0.5. For completeness, the Poisson-Boltzmann prediction 
is provided (dotted line in the upper part of the graph). 



Three values of the coupling constant were considered: 
S — 100, 10 and 0.5. The plots obtained from the equa- 
tion of state (|79|) are shown by dashed lines and the MC 
data 22] are represented by symbols. The accuracy of 
the WSC method is good, surprisingly also for small val- 
ues of S = 10 and 0.5, where the approach is not sup- 
posed to hold. As concerns the (leading term plus the 
first correction) VSC equation of state |23|, correspond- 
ing to our Eq. ^ with^O = S, the plots for 5 = 10 
and 100 are close to the d axis, and far from the Monte 
Carlo data; we consequently do not present them in the 
figure. For S = 0.5, the VSC prediction is in good agree- 



B. Asymmetrically charged plates 

The sequence of ground states for asymmetric like- 
charged plates (0 < C < 1) may be even more complex 
than the one for the symmetric ^ = 1 case; in depen- 
dence on the distance d, the bilayer Wigner crystal can 
involve commensurate as well as incommensurate struc- 
tures of counter-ions. In addition, related work in spher- 
ical geometry [1, has shown that the ground state in 
general breaks local neutrality (the two partners acquire 
an electrical charge, necessarily opposite). The possibil- 
ity of, in principle, an infinite number of irregular struc- 
tures might complicate numerical calculations; we are not 
aware about a work dealing with this subject. 

Fortunately, the same simplification as for the equiva- 
lently charged plates arises at small separations between 
the plates d/b ^ 1, where the lateral lattice spacing b of 
the single Wigner crystal is now given by the requirement 
of the global electroneutrality, as follows: 



ai + a2 



V3 
2 



(80) 



Since the z-coordinates of particles between the plates 
are much smaller than b, we can use the harmonic z- 
expansion of the interaction energy of type (|7ip . where 
only the (irrelevant) constant term reflects the formation 
of some nontrivial bilayer structure. Our task is to derive 
the particle density profile for the energy change from the 
ground state of the form 
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where k=1 — C=l — U2I <y\ and 



N 

E 



4 ^^JR,.-R,|3 

We use the cumulant technique with the one-body Boltz- 
mann factor exp(— The final result for the density 
profile reads 




where 



tl{K) 

<2(k) 



/J* dzze"'*^ 1 



(83) 



(84) 



2d 



(85) 



For example, the density profile p for C = 0.5, S = 86 
and d — 2.68 is depicted in Fig. [3] The dashed curve 
corresponds to the leading SC profile 



po(i) = (l-0- 



1 - e 



-Kd 



(86) 



which is the same in both VSC and WSC theories. For 
the parameters of Fig. |9l the leading order profile reads 



^ ,„ 3 e 
Poiz) = 7 



-z/2 



4 1 



,-1.34 ■ 



(87) 



The WSC profile ((83)) . involving also the first SC correc- 
tion, is represented by the solid curve. The filled circles 
are the MC data of Ref. [l^j. The ratio p/po, which is 
trivially equal to 1 in the leading SC order, is presented 
in the inset of the figure; we see that the first correction 
improves substantially the agreement with MC data. A 
similar conclusion is reached in the case where one plate 
is uncharged (C = 0), see Fig. [TUl for the highest cou- 
pling investigated numerically in Ref. [IJl (S — 86), the 
agreement between the WSC approach and Monte Carlo 
data for the density profile is excellent, and subtle devi- 
ations from the leading order term po are fully captured. 
It can be seen in the inset of Fig. [10] that the agreement 
is no longer quantitative when the coupling parameter 
is decreased by a factor of 10. As may have been an- 
ticipated, the density profile close to the highly charged 



plate located at z = is well accounted for by our treat- 
ment, while the agreement with MC deteriorates when 
approaching the uncharged plate located at z — d. We 
may anticipate that the WSC approach would fare better 
against Monte Carlo at smaller inter-plate separations. 
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FIG. 9. The density profile p for C = 0.5, H = 86 and d = 2.68. 
The dashed curve corresponds to the leading SC profile po 
(|86p . the solid curve also involves the first correction in (|83|l . 
MC data (filled circles) come from Ref. [IJ. The inset shows 
the ratio p/po, for a finer test of the correction to leading 
order po- 
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FIG. 10. Same as in the inset of Fig. [9l for C = 0, and two 
different values of the coupling parameter H. The two plates 
are located at 2 = and z — 2.68. Here, C = means that 
the plate at z — 2.68 is uncharged. The symbols are for the 
Monte Carlo data of Ref. 24]. 



Either of the contact-value relations (|68l) implies the 
same pressure: 

where 

Po - -^(1 + C') + ^(1 - C') coth (^d" 



(89) 
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is the leading SC contribution, already obtained within 
the VSC method in [H, and 



Pi = 



33/4(1 + ^)5/2^^ 



^ 1 coth I ^—^d 



1 



2 / V 2 
is the coefficient of the first l/\/S correction. 
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FIG. 12. Rescaled pressure versus the plate distance for likely 
charged plates with the asymmetry parameter ( = 0.5: The 
dashed curve corresponds to the leading term of the VSC 
theory, which is equivalent to the WSC one (|89[) . The small-d 
expansion of the WSC pressure (|93p is represented by solid 
curves. Filled symbols represent the MC data [SJ] for the 
couplings S = 86 (squares), H = 8.6 (diamonds) and H = 0.32 
(circles in the inset). In the inset, which is a zoom on the 
small distance region, the mean-field Poisson-Boltzmann (PB) 
prediction is also displayed. 



FIG. 11. Oppositely charged plates: The phase boundary 
where P = 0, which discriminates the attractive regime (at 
large distances) from the repulsive one (at small distances). 
The MC data for ( = —0.5 (filled squares) come from Ref. 




While the first correction to the pressure Pi vanishes 
in both limits d — ?> and d — )> cxd, Pq is in general nonzero 
and therefore dominates in these asymptotic regions. Let 
us first consider the large-d limit: 

lim P = Jim Pq = -C^ (91) 

Such a result is correct for oppositely charged plates 
— 1 < C < 0. In that case indeed, for sufficiently dis- 
tant plates, all counter-ions stay in the neighborhood 
of plate 1 and compensate partially its surface charge, 
that is reduced from the bare value aie to \a2\e. We are 
left with a capacitor of opposite surface charges ±1726 
whose dimensionless pressure is attractive and just equal 
to —(^^. In other words, again for large distances, the 
negative counter-ions are expelled from the vicinity of 
the negatively charged plate 2, with a resulting vanish- 
ing charge density p{d). From the contact theorem, this 
implies that the pressure reads P = — C^- Hence, the 
leading SC order (common to VSC and WSC), a pri- 
ori valid at short distances, yields the correct result at 
large distances also. This points to the adequacy of the 
WSC resuh (HHl-jSOl) in the whole range of d values for 
oppositely charged plates, which is consistent with our 
previous analysis about the simple nature of the ground 



state (independent on the inter-plate distance, at vari- 
ance with the C > case). In addition, we emphasize 
that the effect of the first correction coefficient (PH)) is 
very weak. This fact is documented in Fig. [TTJ Each 
solid curve with a fixed asymmetry parameter C < rep- 
resents a phase boundary between the anomalous repul- 
sion of oppositely charged plates at small distances and 
their "natural" attraction at large distances. At S — > 00, 
using the condition Pq = in (15^1) implies the phase 
boundary at [13] 

d*=-2^, S^^ (-1<C<1). (92) 

Considering also the first correction (|M| in ([55]) we see 
in Fig. [TT] that the phase boundary P = is almost 
independent of S, except for very small negative values 
of (. Consequently, the first correction to the leading SC 
behaviour is generically negligible for oppositely charged 
plates. 

On the other hand, the asymptotic result (|9T|) is ap- 
parently physically irrelevant for like-charged plates (0 < 
C < 1). For sufficiently large distances d, the counter- 
ions stay in the neighborhood of both plates 1 and 2 
and a priori neutralize their surface charges, so that the 
asymptotic pressure should vanish. Therefore, for C > Oi 
we cannot expect the same bonus as for ^ < 0, and our 
WSC results dMH-dSOl) hold provided that d < %/S as 
was already the case for C = 1- In addition, the small-d 
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FIG. 13. Phase diagram for like-charged plates with asym- 
metry parameter ( = 0.5. T he p hase boundary given by the 
leading VSC and WSC order [2j is represented by the dashed 
Une. The phase boundary following from our WSC result (|93p 
and (|94|l is represented by the solid curve; for comparison, the 
filled squares are MC data from Ref. [2J] . 
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FIG. 14. The WSC phase boundaries for like-charged plates, 
in the (H, d) plane and for various values of the asymmetry 
parameter 



expansion of the pressure reads 



P = 



i + C , 1 + C 



(1-0^(1 + 



12 



3^(0 



d + 0{(f 



where 



0(0 



(47r)3/2 1 



33/4 Csil + C)^/^ 



(93) 



(94) 



As it should, this is the generalization of the special C = 1 
result ((78l) to all positive asymmetries. 

The plot of the rescaled pressure versus the plate dis- 
tance for likely charged plates with the asymmetry pa- 
rameter C = 0.5 is presented in Fig. [TH The dashed curve 



corresponds to the leading term of the VSC theory, which 
is equivalent to the leading WSC one The small- 

er expansion of the WSC pressure is represented by 
solid curves. The comparison with filled symbols of the 
MC data [l^l shows a good agreement for the coupling 
constants 2 = 86 (squares), S = 8.6 (diamonds) and 
even for relatively small S = 0.32 (circles in the inset). 
The agreement goes somewhat beyond the expected dis- 
tance range of the validity of the expansion ([M)) . but is 
restricted to the small d range. 

The phase diagram for ( = 0.5 is pictured in Fig. [131 
The phase boundary given by the leading 'E, ^ 00 or- 
der of the VSC method [23] is represented by the dashed 
line. As repeatedly emphasized above, it corresponds to 
the leading WSC order as well. The phase boundary fol- 
lowing from our leading plus first correction WSC result 
([M]) and (IMl) is represented by the solid curve; the agree- 
ment with MC data of Ref. [24] (filled squares) is very 
good. The phase boundaries for like-charged plates with 
various values of the asymmetry parameter ^, following 
from our WSC result and (IMl) . are drawn in the 
(E!,d) plane in Fig. [TH It is seen that by decreasing ( 
the anomalous attraction region becomes smaller. 



100 



■T3 50 







1 ' 1 / 




— S = 1000 






■■■■ S = 500 






-- S = 100 






-■■ S = 50 




repulsion 




attraction 


1 1 ■- 1 ; 



0.4 



0.6 



0.8 



FIG. 15. The WSC phase boundaries for like-charged plates, 
in the (^, d) plane and for various values of the coupling con- 
stant H. 



The WSC phase boundaries for like-charged plates, in 
the (C, d) plane and for various values of the coupling 
constant 5, are drawn in Fig. [151 For small values of the 
asymmetry parameter C, e.g. below C, ^ 0.29 for S — 10^, 
we see that the attractive "pocket" disappears. This phe- 
nomenon is entirely driven by the first correction, as in 
revealed by Fig. I16[ which further shows the phase dia- 
gram in the whole range of the asymmetry parameter ^, 
for the coupling constant 5 = 10'^. For comparison, the 
phase boundaries between the repulsion and attractive 
regions in the leading SC order, given by (|92l) . are pic- 
tured by dashed curves. With the corresponding leading 
contribution to the pressure, the attractive region always 
exists. 
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FIG. 16. The WSC phase diagram (sohd curves) in the whole 
range of the asymmetry parameter (, for the coupling con- 
stant H — 10^. For comparison, the phase diagram in the 
leading SC order (|92|) is represented by dashed curves; for 
oppositely charged plates — 1 < C ^ Oi the difference between 
the solid and dashed curves is invisible, due to the already 
pointed out smallness of the first correction for ( < 0. 



IV. CONCLUSION 

In this paper, we have established the mathematical 
grounds for the Wigner Strong CoupUng (WSC) theory 
which describes the strong-coupling regime of counter- 
ions at charged interfaces, starting from the Wigner 
structure formed at zero temperature. The results for 
both likely and oppositely charged plates are in excel- 
lent agreement with Monte Carlo data, which represents 
an improvement over the previously proposed Virial SC 
approach. By construction, our expansion should be 
more reliable the larger the coupling parameter S, but 
we found that it remains trustworthy for intermediate 
values of the coupling constant (say S = 100), and in 
some cases down to S = 10 or 20. 

The geometries studied are those of one or two planar 
interfaces. An important remark is that the leading re- 
sults in the SC expansion follow from a single counter-ion 
picture because the dominant (linear) electric potential 
stems from the plate only; the contribution due to the 
interaction with other counter-ions on the same plate is 
harmonic and therefore sub-dominant. As a consequence, 
the leading terms of the VSC and WSC theories coincide. 
This fact has been outlined on several occasions, but can 



nevertheless not be considered as a general statement. 
Indeed, the situation changes for a curved (say, cylin- 
drical or spherical) wall surfaces since then the interac- 
tions of an ion with other counter-ions contribute to the 
dominant field, no matter how close to the interface this 
ion can be. This is why the leading ion profile around a 
charged cylinder or sphere will in general differ from that 
obtained within the original VSC approach [27]. Inclu- 
sion of curvature effects in the WSC treatment is a task 
for the future. In the present work, we have also assumed 
that the charges on the plates are uniformly smeared, 
which opens the way to the powerful use of the contact 
theorem to obtain the pressure. As a consequence, the 
interesting case of discrete fixed charges on the plates 
[45l - l48l |. is beyond the scope of the present analysis. 

A generalization of the formalism to quantum statis- 
tical systems of counter-ions is straightforward: Vibra- 
tions of counter-ions around their Wigner-lattice posi- 
tions possess energy spectrum of quantized harmonic os- 
cillators. Another perspective is to formulate a strong- 
coupling theory valid for an arbitrary distance between 
the plates. Indeed, both the original Virial SC and the 
present Wigner SC theories are so far limited, in the two 
plate case, to the regime d <€i 2^/^, which means that 
the inter plate distance should be smaller than the lat- 
tice spacing a in the underlying Wigner crystal (up to 
an irrelevant prefactor, the quantities a and b introduced 
in this article refer to the same length). It is impor- 
tant to emphasize here that the limitation d <C ^^^^ is 
not intrinsic to the strong coupling limit, but is a tech- 
nical requirement that should be enforced to allow for 
the validity of the single particle picture, and subsequent 
higher order corrections as worked out here. Perform- 
ing the SC expansion for distances d ^ S^/^ requires to 
bypass the single particle picture, which is a challenging 
goal. Finally, in view of possible applications to real col- 
loidal systems, it seems important to account for the low 
dielectric constant of colloidal particles, taking due ac- 
count of image charge effects ^25i . i49i] . Work along these 
lines is in progress. 
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